#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>

#include "MSTK.h"


#ifdef __cplusplus
extern "C" {
#endif

int MESH_Vol_AddGhost(Mesh_ptr mesh, Mesh_ptr submesh, int ring) {
  int iv, i, j, idx, vertex_id, edge_id, face_id;
  int nv, ne, nf, nvf, nfv, nfe;
  MVertex_ptr mv_new, *fverts, global_mv, local_mv, global_mfvertices;
  MEdge_ptr local_me, global_me, local_mfedges;
  MFace_ptr local_mf, global_mf, mf_new;
  List_ptr global_mvfaces, global_mfedges;

  /* basic mesh information */
  nv = MESH_Num_Vertices(submesh);
  ne = MESH_Num_Edges(submesh);
  nf = MESH_Num_Faces(submesh);

  /* Mark the list of global entities in this submesh*/
  List_ptr par_vertices=List_New(nv);
  int mkvid = MSTK_GetMarker();
  idx = 0;
  while((local_mv = MESH_Next_Vertex(submesh,&idx))) {
    vertex_id = MEnt_ID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    List_Add(par_vertices,global_mv);
    MEnt_Mark(global_mv,mkvid);
  }
  fprintf(stdout,"number of vertices: %d, marked: %d\n",nv,List_Num_Entries(par_vertices));
  List_ptr par_edges = List_New(ne);
  int mkeid = MSTK_GetMarker();
  idx = 0;
  while((local_me = MESH_Next_Edge(submesh,&idx))) {
    edge_id = MEnt_ID(local_me);
    global_me = MESH_EdgeFromID(mesh,edge_id);
    List_Add(par_edges,global_me);
    MEnt_Mark(global_me,mkeid);
  }
  fprintf(stdout,"number of edges: %d, marked: %d\n",ne,List_Num_Entries(par_edges));
  List_ptr par_faces=List_New(nf);
  int mkfid = MSTK_GetMarker();  
  idx = 0;
  while((local_mf = MESH_Next_Face(submesh,&idx))) {
    face_id = MEnt_ID(local_mf);
    global_mf = MESH_FaceFromID(mesh,face_id);
    List_Add(par_faces,global_mf);
    MEnt_Mark(global_mf,mkfid);
  }
  fprintf(stdout,"number of faces: %d, marked: %d\n",nf,List_Num_Entries(par_faces));


  /*Loop through sub vertices, and see its neighboring faces */
  int index_face = 0;
  for(iv = 0; iv < nv; iv++) {
    /* get global face id and loop through global mesh
     here, can not use List par_vertices, it will be dynamically added */
    local_mv = MESH_Vertex(submesh,iv);
    /* only check vertices on the boundary */
    if(MV_PType(local_mv) != PBOUNDARY)
      continue;
    
    vertex_id = MEnt_ID(local_mv);
    global_mv = MESH_VertexFromID(mesh,vertex_id);
    /* get global neighbor faces of this vertex */
    global_mvfaces = MV_Faces(global_mv);
    nvf = List_Num_Entries(global_mvfaces);
    for(i = 0; i < nvf; i++) {
      global_mf = List_Entry(global_mvfaces,i);
      /* if it is not marked, add this face into submesh */
      if (!MEnt_IsMarked(global_mf,mkfid)) {
	/* loop through the global vertices of this face */
	global_mfvertices = MF_Vertices(global_mf,1,0);
	nfv = List_Num_Entries(global_mfvertices);
	fverts = (MVertex_ptr *)MSTK_malloc((nfv)*sizeof(MVertex_ptr));
	for(j = 0; j < nfv; j++) {
	  global_mv = List_Entry(global_mfvertices,j);
	    /* if this vertex is already marked, use its local counterpart */
	  if(MEnt_IsMarked(global_mv,mkvid)) {
	    fverts[j] = MESH_VertexFromID(submesh,MV_ID(global_mv));
	    if(MV_PType(local_mv) != PGHOST)
	      MV_Set_PType(local_mv,PBOUNDARY);
	  }
	  /* otherwise, add this vertex into submesh */
	  else {
	    List_Add(par_vertices,global_mv);
	    MEnt_Mark(global_mv,mkvid);
	    mv_new = MV_New(submesh);
	    PMV_Duplicate(global_mv,mv_new);
	    /* ptype is not copied in PMV_Copy */
	    MV_Set_PType(mv_new,PGHOST);
	    fverts[j] = mv_new;
	  }
	}
	/* add face */
	List_Add(par_faces,global_mf);
	index_face++;
	MEnt_Mark(global_mf,mkfid);
	mf_new = MF_New(submesh);
	MF_Set_Vertices(mf_new,nfv,fverts);
	MF_Set_PType(mf_new,PGHOST);
	PMF_Duplicate(global_mf,mf_new);
	if (fverts) MSTK_free(fverts);
	List_Delete(global_mfvertices);
	/* add edges */
	local_mfedges = MF_Edges(mf_new,1,0);
	/* after call MF_Set_Vertices(), newly added edges ID were reset 
	 here we recover the global ID */
	global_mfedges = MF_Edges(global_mf,1,0);
	nfe = List_Num_Entries(local_mfedges);
	for (j = 0; j < nfe; j++) {
	  local_me = List_Entry(local_mfedges,j);
	  global_me = List_Entry(global_mfedges,j);
	  ME_Set_ID(local_me,MV_ID(global_me));
	  /* if not marked, mark as ghost edge */
	  if(!MEnt_IsMarked(global_me, mkeid)) {
	    List_Add(par_edges,global_me);
	    MEnt_Mark(global_me,mkeid);
	    ME_Set_PType(local_me,PGHOST);
	    PME_Duplicate(global_me,local_me);
	  }
	}
	List_Delete(local_mfedges);
	List_Delete(global_mfedges);
      }
    }
    List_Delete(global_mvfaces);
  }
  fprintf(stdout,"num of faces added: %d\n",index_face);
  fprintf(stdout,"num of faces on submesh: %d\n",MESH_Num_Faces(submesh));
  List_Unmark(par_vertices,mkvid);
  List_Unmark(par_edges,mkeid);
  List_Unmark(par_faces,mkfid);
  MSTK_FreeMarker(mkvid);
  MSTK_FreeMarker(mkeid);
  MSTK_FreeMarker(mkfid);
  fprintf(stdout,"number of ghost vertices: %d\n",List_Num_Entries(par_vertices));
  fprintf(stdout,"number of ghost edges: %d\n",List_Num_Entries(par_edges));
  fprintf(stdout,"number of ghost faces: %d\n",List_Num_Entries(par_faces));
  List_Delete(par_vertices);
  List_Delete(par_edges);
  List_Delete(par_faces);
  return 1;
}
  
#ifdef __cplusplus
}
#endif

